#Pathway analysis for Proteomic data using CePa
setwd("/Users/emmatimminsschiffman/Documents/Dissertation/proteomics/DB post-genome/OT expression comparison/pathway analysis")

#First will do gene set analysis (GSA) method by using an expression matrix and phenotype label
#create sample label for pCO2 analysis (it seems that the function can only handle one treatment and one control so will have to divide up treatments)
treatments<-sampleLabel(c(rep("HighCO2", 8), rep("LowCO2", 8)), treatment="HighCO2", control="LowCO2")

#expression matrix of all proteins (proteins are grouped by SPID)
exp.mat<-read.csv('all protein spid expression for CePa.csv', header=T, row.names=1)

#create pathway list
data(PID.db)

#run cepa for GSA
pathway.gsa<-cepa.all(mat=exp.mat, label=treatments, pc=PID.db$NCI)
report(pathway.gsa, adj.method='BH', cutoff=0.1, only.sig=TRUE)
#no significant pathways

#Effect of MS across both pCO2
treatment.ms<-sampleLabel(c(rep("NS", 4), rep("MS", 4), rep("NS", 4), rep("MS", 4)), treatment="MS", control="NS")

pathway.gsa2<-cepa.all(mat=exp.mat, label=treatment.ms, pc=PID.db$NCI)
plot(pathway.gsa2)

#effect of MS only at ambient pCO2
treatment.ms2<-sampleLabel(c(rep("NS",4), rep("MS", 4)), treatment="MS", control="NS")

pathway.gsa3<-cepa.all(mat=exp.mat[9:16], label=treatment.ms2, pc=PID.db$NCI)
plot(pathway.gsa3)

#Second will do Over-representation analysis (ORA)
setwd("/Users/emmatimminsschiffman/Documents/Dissertation/proteomics/DB post-genome/OT expression comparison")
bkgrnd<-read.table('DAVID background at least 10 specs.txt')
protein.list<-read.table('DAVID_input_with_SPID.txt', sep="\t", header=TRUE)

pathway.ora<-cepa.all(dif=protein.list$High, bk=bkgrnd, pc=PID.db$NCI)
report(pathway.ora, adj.method='BH', cutoff=0.1)
#no.sig

pathway.ora2<-cepa.all(dif=protein.list$High.MS, bk=bkgrnd, pc=PID.db$NCI)
plot(pathway.ora2)

pathway.ora3<-cepa.all(dif=protein.list$Ambient, bk=bkgrnd, pc=PID.db$NCI)
plot(pathway.ora3)

pathway.ora4<-cepa.all(dif=protein.list$Ambient.MS, bk=bkgrnd, pc=PID.db$NCI)
plot(pathway.ora4)